# install.packages("MatchIt")
# install.packages("readstata13")
# install.packages("MASS")
# install.packages("dplyr")
# install.packages("survey")
# install.packages("tidyr")
# install.packages("ggplot2")
# install.packages("haven")
# install.packages("stats")
# install.packages("Matrix")
# install.packages("data.table")

library(MatchIt)
library(readstata13)
library(MASS)
library(dplyr)
library(survey)
library(tidyr)
library(ggplot2)
library(haven)
library(stats)
library(Matrix)
library(data.table)

# Data
data <- read_dta("../../data/processed/ces/sensitivity_checks/add_housing/rents_hh_2022_prematch.dta")

# Keep only for positive values of rent
data <- data %>%
  filter(rent_w1 > 0)

#Keep variables to input in the model
data_model = subset(data, 
                    select = c(l_hhsize, wall_material, roof_material,
                               floor_material, drinking_water, latrine_type,
                               lighting, mobile, internet, cooking, 
                               log_time_to_water, caste, television, laptop,
                               air_conditioner, car, refrigerator,
                               washing_machine, nic_agro, nic_industry, 
                               nic_services, nic_public, hindu, l_rent, pwt,
                               state, sector, sss, survey_month, renter, hhid))

rm(data)
#data_model$hhidseq=1:nrow(data_model)
# Define categorical variables
data_model <- data_model %>%
  mutate(hhid=as.numeric(hhid),
         state = as.factor(state),
         wall_material = as.factor(wall_material),
         roof_material = as.factor(roof_material),
         floor_material = as.factor(floor_material),
         latrine_type = as.factor(latrine_type),
         drinking_water = as.factor(drinking_water),
         lighting = as.factor(lighting),
         cooking = as.factor(cooking),
         caste = as.factor(caste))

#Remove NAs
data_model=na.omit(data_model)


#Full model for variable selection
full_model <- glm(renter~ . -state -sector -sss -l_rent -pwt -survey_month -hhid, 
               data = data_model, family = binomial)


#Variable selection via AIC based stepwise
sel_model_c <- stepAIC(full_model, direction = "both")

#Model summary
summary(sel_model_c)

# Get the list of selected variables
sel_variables_c <- names(coef(sel_model_c))[-1]  # Exclude the intercept

#data_model2 <- data_model %>% 
data_matrix = model.matrix(~ . -state -sector -sss -l_rent -pwt -survey_month, 
                          data = data_model)

data_matrix <- data_matrix[, 
          colnames(data_matrix) %in% c(sel_variables_c,"hhid")]

data_df=as.data.frame(as.matrix(data_matrix))

rm(data_matrix)

data_df=merge(data_df,data_model[,
                c("hhid","state","sector","renter")],by="hhid")

# Add division and urban/rural for exact matching
sel_variables_c = c(sel_variables_c,c("state", "sector"))

# Create the formula for matching using the selected variables
matching_formula_c <- as.formula(paste("renter ~", 
                                       paste(sel_variables_c,  collapse = " + ")))

# Exact matching on division and urban rural
matchit_out_c <- matchit(
  matching_formula_c,
  data = data_df,
  distance = "glm",
  method = "nearest",
  exact = ~ state + sector,
  ratio = NULL,  # Allows for repetition (many-to-one matching)
  replace = TRUE,  # Allow treated units to be matched with multiple controls
  discard="both"
)

# Summary of the matching result
summary(matchit_out_c)
matchit_out_c$model

# Extract the matched dataset
matched_data_c <- match.data(matchit_out_c)

matched_data_c$hhid = as.character(matched_data_c$hhid)

# Save dataset
write_dta(matched_data_c, "../../data/processed/ces/sensitivity_checks/add_housing/train_R.dta")
